#include <cstdlib>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <ctime>
#include <cmath>
#include <cstring>

using namespace std;

#include "table_io.h"

//****************************************************************************80

char ch_cap(char ch)

//****************************************************************************80
//
//  Purpose:
//
//    CH_CAP capitalizes a single character.
//
//  Discussion:
//
//    This routine should be equivalent to the library "toupper" function.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    19 July 1998
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, char CH, the character to capitalize.
//
//    Output, char CH_CAP, the capitalized character.
//
{
    if (97 <= ch && ch <= 122) {
        ch = ch - 32;
    }

    return ch;
}
//****************************************************************************80

bool ch_eqi(char ch1, char ch2)

//****************************************************************************80
//
//  Purpose:
//
//    CH_EQI is true if two characters are equal, disregarding case.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    13 June 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, char CH1, CH2, the characters to compare.
//
//    Output, bool CH_EQI, is true if the two characters are equal,
//    disregarding case.
//
{
    if (97 <= ch1 && ch1 <= 122) {
        ch1 = ch1 - 32;
    }
    if (97 <= ch2 && ch2 <= 122) {
        ch2 = ch2 - 32;
    }

    return ( ch1 == ch2);
}
//****************************************************************************80

int ch_to_digit(char ch)

//****************************************************************************80
//
//  Purpose:
//
//    CH_TO_DIGIT returns the integer value of a base 10 digit.
//
//  Example:
//
//     CH  DIGIT
//    ---  -----
//    '0'    0
//    '1'    1
//    ...  ...
//    '9'    9
//    ' '    0
//    'X'   -1
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    13 June 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, char CH, the decimal digit, '0' through '9' or blank are legal.
//
//    Output, int CH_TO_DIGIT, the corresponding integer value.  If the 
//    character was 'illegal', then DIGIT is -1.
//
{
    int digit;

    if ('0' <= ch && ch <= '9') {
        digit = ch - '0';
    } else if (ch == ' ') {
        digit = 0;
    } else {
        digit = -1;
    }

    return digit;
}
//****************************************************************************80

int file_column_count(string filename)

//****************************************************************************80
//
//  Purpose:
//
//    FILE_COLUMN_COUNT counts the columns in the first line of a file.
//
//  Discussion:
//
//    The file is assumed to be a simple text file.
//
//    Most lines of the file are presumed to consist of COLUMN_NUM words, 
//    separated by spaces.  There may also be some blank lines, and some 
//    comment lines, which have a "#" in column 1.
//
//    The routine tries to find the first non-comment non-blank line and
//    counts the number of words in that line.
//
//    If all lines are blanks or comments, it goes back and tries to analyze
//    a comment line.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    05 July 2009
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, string FILENAME, the name of the file.
//
//    Output, int FILE_COLUMN_COUNT, the number of columns assumed 
//    to be in the file.
//
{
    int column_num;
    ifstream input;
    bool got_one;
    string text;
    //
    //  Open the file.
    //
    input.open(filename.c_str());

    if (!input) {
        column_num = -1;
        cerr << "\n";
        cerr << "FILE_COLUMN_COUNT - Fatal error!\n";
        cerr << "  Could not open the file:\n";
        cerr << "  \"" << filename << "\"\n";
        return column_num;
    }
    //
    //  Read one line, but skip blank lines and comment lines.
    //
    got_one = false;

    for (;;) {
        getline(input, text);

        if (input.eof()) {
            break;
        }

        if (s_len_trim(text) <= 0) {
            continue;
        }

        if (text[0] == '#' || text[0] == '/' || text[0] == '!') {
            continue;
        }
        got_one = true;
        break;
    }

    if (!got_one) {
        input.close();

        input.open(filename.c_str());

        for (;;) {
            input >> text;

            if (input.eof()) {
                break;
            }

            if (s_len_trim(text) == 0) {
                continue;
            }
            got_one = true;
            break;
        }
    }

    input.close();

    if (!got_one) {
        cerr << "\n";
        cerr << "FILE_COLUMN_COUNT - Warning!\n";
        cerr << "  The file does not seem to contain any data.\n";
        return -1;
    }

    column_num = s_word_count(text);

    return column_num;
}
//****************************************************************************80

int file_row_count(string input_filename)

//****************************************************************************80
//
//  Purpose:
//
//    FILE_ROW_COUNT counts the number of row records in a file.
//
//  Discussion:
//
//    It does not count lines that are blank, or that begin with a
//    comment symbol '#'.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 February 2009
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, string INPUT_FILENAME, the name of the input file.
//
//    Output, int FILE_ROW_COUNT, the number of rows found.
//
{
    int comment_num;
    ifstream input;
    string line;
    int record_num;
    int row_num;

    row_num = 0;
    comment_num = 0;
    record_num = 0;

    input.open(input_filename.c_str());

    if (!input) {
        cerr << "\n";
        cerr << "FILE_ROW_COUNT - Fatal error!\n";
        cerr << "  Could not open the input file: \"" << input_filename << "\"\n";
        return (-1);
    }

    for (;;) {
        getline(input, line);

        if (input.eof()) {
            break;
        }

        record_num = record_num + 1;

        if (line[0] == '#' || line[0] == '/' || line[0] == '!') {
            comment_num = comment_num + 1;
            continue;
        }

        if (s_len_trim(line) == 0) {
            comment_num = comment_num + 1;
            continue;
        }

        row_num = row_num + 1;

    }

    input.close();

    return row_num;
}
//****************************************************************************80

int i4_log_10(int i)

//****************************************************************************80
//
//  Purpose:
//
//    I4_LOG_10 returns the integer part of the logarithm base 10 of ABS(X).
//
//  Example:
//
//        I  I4_LOG_10
//    -----  --------
//        0    0
//        1    0
//        2    0
//        9    0
//       10    1
//       11    1
//       99    1
//      100    2
//      101    2
//      999    2
//     1000    3
//     1001    3
//     9999    3
//    10000    4
//
//  Discussion:
//
//    I4_LOG_10 ( I ) + 1 is the number of decimal digits in I.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    04 January 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int I, the number whose logarithm base 10 is desired.
//
//    Output, int I4_LOG_10, the integer part of the logarithm base 10 of
//    the absolute value of X.
//
{
    int i_abs;
    int ten_pow;
    int value;

    if (i == 0) {
        value = 0;
    } else {
        value = 0;
        ten_pow = 10;

        i_abs = abs(i);

        while (ten_pow <= i_abs) {
            value = value + 1;
            ten_pow = ten_pow * 10;
        }

    }

    return value;
}
//****************************************************************************80

int i4_max(int i1, int i2)

//****************************************************************************80
//
//  Purpose:
//
//    I4_MAX returns the maximum of two I4's.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    13 October 1998
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int I1, I2, are two integers to be compared.
//
//    Output, int I4_MAX, the larger of I1 and I2.
//
{
    int value;

    if (i2 < i1) {
        value = i1;
    } else {
        value = i2;
    }
    return value;
}
//****************************************************************************80

int i4_min(int i1, int i2)

//****************************************************************************80
//
//  Purpose:
//
//    I4_MIN returns the minimum of two I4's.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    13 October 1998
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int I1, I2, two integers to be compared.
//
//    Output, int I4_MIN, the smaller of I1 and I2.
//
{
    int value;

    if (i1 < i2) {
        value = i1;
    } else {
        value = i2;
    }
    return value;
}
//****************************************************************************80

int *i4mat_border_add(int m, int n, int table[])

//****************************************************************************80
//
//  Purpose:
//
//    I4MAT_BORDER_ADD adds a "border" to an I4MAT.
//
//  Discussion:
//
//    An I4MAT is an array of I4's.
//
//    We suppose the input data gives values of a quantity on nodes
//    in the interior of a 2D grid, and we wish to create a new table
//    with additional positions for the nodes that would be on the
//    border of the 2D grid.
//
//                  0 0 0 0 0 0
//      * * * *     0 * * * * 0
//      * * * * --> 0 * * * * 0
//      * * * *     0 * * * * 0
//                  0 0 0 0 0 0
//
//    The illustration suggests the situation in which a 3 by 4 array
//    is input, and a 5 by 6 array is to be output.
//
//    The old data is shifted to its correct positions in the new array.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    25 January 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int M, the spatial dimension.
//
//    Input, int N, the number of points.
//
//    Input, int TABLE[M*N], the table data.
//
//    Output, int TABLE2[(M+2)*(N+2)], the augmented table data.
//
{
    int i;
    int j;
    int *table2;

    table2 = new int[(m + 2)*(n + 2)];

    for (j = 0; j < n + 2; j++) {
        for (i = 0; i < m + 2; i++) {
            if (i == 0 || i == m + 1 || j == 0 || j == n + 1) {
                table2[i + j * (m + 2)] = 0;
            } else {
                table2[i + j * (m + 2)] = table[(i - 1)+(j - 1) * m];
            }
        }
    }

    return table2;
}
//****************************************************************************80

int *i4mat_border_cut(int m, int n, int table[])

//****************************************************************************80
//
//  Purpose:
//
//    I4MAT_BORDER_CUT cuts the "border" of an I4MAT.
//
//  Discussion:
//
//    An I4MAT is an array of I4's.
//
//    We suppose the input data gives values of a quantity on nodes
//    on a 2D grid, and we wish to create a new table corresponding only
//    to those nodes in the interior of the 2D grid.
//
//      0 0 0 0 0 0
//      0 * * * * 0    * * * *
//      0 * * * * 0 -> * * * *
//      0 * * * * 0    * * * *
//      0 0 0 0 0 0
//
//    The illustration suggests the situation in which a 5 by 6 array
//    is input, and a 3 by 4 array is to be output.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    25 January 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int M, the spatial dimension.
//
//    Input, int N, the number of points.
//
//    Input, int TABLE[M*N], the table data.
//
//    Output, int TABLE2[(M-2)*(N-2)], the "interior" table data.
//
{
    int i;
    int j;
    int *table2;

    if (m <= 2 || n <= 2) {
        return NULL;
    }

    table2 = new int[(m - 2)*(n - 2)];

    for (j = 0; j < n - 2; j++) {
        for (i = 0; i < m - 2; i++) {
            table2[i + j * (m - 2)] = table[(i + 1)+(j + 1) * m];
        }
    }

    return table2;
}
//****************************************************************************80

int *i4mat_data_read(string input_filename, int m, int n)

//****************************************************************************80
//
//  Purpose:
//
//    I4MAT_DATA_READ reads data from an I4MAT file.
//
//  Discussion:
//
//    An I4MAT is an array of I4's.
//
//    The file is assumed to contain one record per line.
//
//    Records beginning with '#' are comments, and are ignored.
//    Blank lines are also ignored.
//
//    Each line that is not ignored is assumed to contain exactly (or at least)
//    M real numbers, representing the coordinates of a point.
//
//    There are assumed to be exactly (or at least) N such records.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 February 2009
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, string INPUT_FILENAME, the name of the input file.
//
//    Input, int M, the number of spatial dimensions.
//
//    Input, int N, the number of points.  The program
//    will stop reading data once N values have been read.
//
//    Output, int I4MAT_DATA_READ[M*N], the table data.
//
{
    bool error;
    ifstream input;
    int i;
    int j;
    string line;
    int *table;
    int *x;

    input.open(input_filename.c_str());

    if (!input) {
        cerr << "\n";
        cerr << "I4MAT_DATA_READ - Fatal error!\n";
        cerr << "  Could not open the input file: \"" << input_filename << "\"\n";
        return NULL;
    }

    table = new int[m * n];

    x = new int[m];

    j = 0;

    while (j < n) {
        getline(input, line);

        if (input.eof()) {
            break;
        }

        if (line[0] == '#' || line[0] == '/' || line[0] == '!' || s_len_trim(line) == 0) {
            continue;
        }

        error = s_to_i4vec(line, m, x);

        if (error) {
            continue;
        }

        for (i = 0; i < m; i++) {
            //      table[i+j*m] = x[i];
            table[i * n + j] = x[i];
        }
        j = j + 1;

    }

    input.close();

    delete [] x;

    return table;
}
//****************************************************************************80

void i4mat_header_read(string input_filename, int *m, int *n)

//****************************************************************************80
//
//  Purpose:
//
//    I4MAT_HEADER_READ reads the header from an I4MAT file.
//
//  Discussion:
//
//    An I4MAT is an array of I4's.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 February 2009
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, string INPUT_FILENAME, the name of the input file.
//
//    Output, int *M, the number of spatial dimensions.
//
//    Output, int *N, the number of points
//
{
    *m = file_column_count(input_filename);

    if (*m <= 0) {
        cerr << "\n";
        cerr << "I4MAT_HEADER_READ - Fatal error!\n";
        cerr << "  FILE_COLUMN_COUNT failed.\n";
        *n = -1;
        return;
    }

    *n = file_row_count(input_filename);

    if (*n <= 0) {
        cerr << "\n";
        cerr << "I4MAT_HEADER_READ - Fatal error!\n";
        cerr << "  FILE_ROW_COUNT failed.\n";
        return;
    }

    return;
}
//****************************************************************************80

int *i4mat_indicator(int m, int n)

//****************************************************************************80
//
//  Purpose:
//
//    I4MAT_INDICATOR sets up an "indicator" I4MAT.
//
//  Discussion:
//
//    An I4MAT is an array of I4's.
//
//    The value of each entry suggests its location, as in:
//
//      11  12  13  14
//      21  22  23  24
//      31  32  33  34
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    25 January 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int M, the number of rows of the matrix.
//    M must be positive.
//
//    Input, int N, the number of columns of the matrix.
//    N must be positive.
//
//    Output, int TABLE[M*N], the indicator matrix.
//
{
    int fac;
    int i;
    int j;
    int *table;

    table = new int[m * n];

    fac = (int) pow(10.0, (i4_log_10(n) + 1));

    for (i = 1; i <= m; i++) {
        for (j = 1; j <= n; j++) {
            table[i - 1 + (j - 1) * m] = fac * i + j;
        }
    }

    return table;
}
//****************************************************************************80

void i4mat_print(int m, int n, int a[], string title)

//****************************************************************************80
//
//  Purpose:
//
//    I4MAT_PRINT prints an I4MAT, with an optional title.
//
//  Discussion:
//
//    An I4MAT is an array of I4's.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    30 April 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int M, the number of rows in A.
//
//    Input, int N, the number of columns in A.
//
//    Input, int A[M*N], the M by N matrix.
//
//    Input, string TITLE, a title to be printed.
//
{

    i4mat_print_some(m, n, a, 1, 1, m, n, title);

    return;
}
//****************************************************************************80

void i4mat_print_some(int m, int n, int a[], int ilo, int jlo, int ihi,
        int jhi, string title)

//****************************************************************************80
//
//  Purpose:
//
//    I4MAT_PRINT_SOME prints some of an I4MAT.
//
//  Discussion:
//
//    An I4MAT is an array of I4's.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    09 April 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int M, the number of rows of the matrix.
//    M must be positive.
//
//    Input, int N, the number of columns of the matrix.
//    N must be positive.
//
//    Input, int A[M*N], the matrix.
//
//    Input, int ILO, JLO, IHI, JHI, designate the first row and
//    column, and the last row and column to be printed.
//
//    Input, string TITLE, a title for the matrix.
{
#define INCX 10

    int i;
    int i2hi;
    int i2lo;
    int j;
    int j2hi;
    int j2lo;

    if (0 < s_len_trim(title)) {
        cout << "\n";
        cout << title << "\n";
    }
    //
    //  Print the columns of the matrix, in strips of INCX.
    //
    for (j2lo = jlo; j2lo <= jhi; j2lo = j2lo + INCX) {
        j2hi = j2lo + INCX - 1;
        j2hi = i4_min(j2hi, n);
        j2hi = i4_min(j2hi, jhi);

        cout << "\n";
        //
        //  For each column J in the current range...
        //
        //  Write the header.
        //
        cout << "  Col: ";
        for (j = j2lo; j <= j2hi; j++) {
            cout << setw(6) << j << "  ";
        }
        cout << "\n";
        cout << "  Row\n";
        cout << "  ---\n";
        //
        //  Determine the range of the rows in this strip.
        //
        i2lo = i4_max(ilo, 1);
        i2hi = i4_min(ihi, m);

        for (i = i2lo; i <= i2hi; i++) {
            //
            //  Print out (up to INCX) entries in row I, that lie in the current strip.
            //
            cout << setw(5) << i << "  ";
            for (j = j2lo; j <= j2hi; j++) {
                cout << setw(6) << a[i - 1 + (j - 1) * m] << "  ";
            }
            cout << "\n";
        }
    }

    return;
#undef INCX
}
//****************************************************************************80

int *i4mat_read(string input_filename, int *m, int *n)

//****************************************************************************80
//
//  Purpose:
//
//    I4MAT_READ reads the information from an I4MAT file.
//
//  Discussion:
//
//    An I4MAT is an array of I4's.
//
//    The file is assumed to contain one record per line.
//
//    Records beginning with '#' are comments, and are ignored.
//    Blank lines are also ignored.
//
//    Each line that is not ignored is assumed to contain exactly (or at least)
//    M real numbers, representing the coordinates of a point.
//
//    There are assumed to be exactly (or at least) N such records.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 February 2009
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, string INPUT_FILENAME, the name of the input file.
//
//    Output, int *M, the number of spatial dimensions.
//
//    Output, int *N, the number of points.  The program
//    will stop reading data once N values have been read.
//
//    Output, int I4MAT_READ[M*N], the table data.
//
{
    int *table;

    i4mat_header_read(input_filename, m, n);

    table = i4mat_data_read(input_filename, *m, *n);

    return table;
}
//****************************************************************************80

void i4mat_write(string output_filename, int m, int n, int table[])

//****************************************************************************80
//
//  Purpose:
//
//    I4MAT_WRITE writes an I4MAT file with no header.
//
//  Discussion:
//
//    An I4MAT is an array of I4's.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    01 June 2009
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, string OUTPUT_FILENAME, the output filename.
//
//    Input, int M, the spatial dimension.
//
//    Input, int N, the number of points.
//
//    Input, int TABLE[M*N], the table data.
//
{
    int i;
    int j;
    ofstream output;
    //
    //  Open the file.
    //
    output.open(output_filename.c_str());

    if (!output) {
        cerr << "\n";
        cerr << "I4MAT_WRITE - Fatal error!\n";
        cerr << "  Could not open the output file.\n";
        return;
    }
    //
    //  Write the data.
    //
    for (j = 0; j < n; j++) {
        for (i = 0; i < m; i++) {
            output << "  " << setw(10) << table[i + j * m];
        }
        output << "\n";
    }
    //
    //  Close the file.
    //
    output.close();

    return;
}
//****************************************************************************80

double *r8mat_data_read(string input_filename, int m, int n)

//****************************************************************************80
//
//  Purpose:
//
//    R8MAT_DATA_READ reads the data from an R8MAT file.
//
//  Discussion:
//
//    An R8MAT is an array of R8's.
//
//    The file is assumed to contain one record per line.
//
//    Records beginning with '#' are comments, and are ignored.
//    Blank lines are also ignored.
//
//    Each line that is not ignored is assumed to contain exactly (or at least)
//    M real numbers, representing the coordinates of a point.
//
//    There are assumed to be exactly (or at least) N such records.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 February 2009
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, string INPUT_FILENAME, the name of the input file.
//
//    Input, int M, the number of spatial dimensions.
//
//    Input, int N, the number of points.  The program
//    will stop reading data once N values have been read.
//
//    Output, double R8MAT_DATA_READ[M*N], the table data.
//
{
    bool error;
    ifstream input;
    int i;
    int j;
    string line;
    double *table;
    double *x;

    input.open(input_filename.c_str());

    if (!input) {
        cerr << "\n";
        cerr << "R8MAT_DATA_READ - Fatal error!\n";
        cerr << "  Could not open the input file: \"" << input_filename << "\"\n";
        return NULL;
    }

    table = new double[m * n];

    x = new double[m];

    j = 0;

    while (j < n) {
        getline(input, line);

        if (input.eof()) {
            break;
        }

        if (line[0] == '#' || line[0] == '/' || line[0] == '!' || s_len_trim(line) == 0) {
            continue;
        }

        error = s_to_r8vec(line, m, x);

        if (error) {
            continue;
        }

        for (i = 0; i < m; i++) {
            //      table[i+j*m] = x[i];
            table[i * n + j] = x[i];
        }
        j = j + 1;

    }

    input.close();

    delete [] x;

    return table;
}
//****************************************************************************80

float *r4mat_data_read(string input_filename, int m, int n)

//****************************************************************************80
//
//  Purpose:
//
//    R8MAT_DATA_READ reads the data from an R8MAT file.
//
//  Discussion:
//
//    An R8MAT is an array of R8's.
//
//    The file is assumed to contain one record per line.
//
//    Records beginning with '#' are comments, and are ignored.
//    Blank lines are also ignored.
//
//    Each line that is not ignored is assumed to contain exactly (or at least)
//    M real numbers, representing the coordinates of a point.
//
//    There are assumed to be exactly (or at least) N such records.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 February 2009
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, string INPUT_FILENAME, the name of the input file.
//
//    Input, int M, the number of spatial dimensions.
//
//    Input, int N, the number of points.  The program
//    will stop reading data once N values have been read.
//
//    Output, double R8MAT_DATA_READ[M*N], the table data.
//
{
    bool error;
    ifstream input;
    int i;
    int j;
    string line;
    float *table;
    double *x;

    input.open(input_filename.c_str());

    if (!input) {
        cerr << "\n";
        cerr << "R8MAT_DATA_READ - Fatal error!\n";
        cerr << "  Could not open the input file: \"" << input_filename << "\"\n";
        return NULL;
    }

    table = new float[m * n];

    x = new double[m];

    j = 0;

    while (j < n) {
        getline(input, line);

        if (input.eof()) {
            break;
        }

        if (line[0] == '#' || line[0] == '/' || line[0] == '!' || s_len_trim(line) == 0) {
            continue;
        }

        error = s_to_r8vec(line, m, x);

        if (error) {
            continue;
        }

        for (i = 0; i < m; i++) {
            //      table[i+j*m] = x[i];
            table[i * n + j] = x[i];
        }
        j = j + 1;

    }

    input.close();

    delete [] x;

    return table;
}
//****************************************************************************80

void r8mat_header_read(string input_filename, int *m, int *n)

//****************************************************************************80
//
//  Purpose:
//
//    R8MAT_HEADER_READ reads the header from an R8MAT file.
//
//  Discussion:
//
//    An R8MAT is an array of R8's.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 February 2009
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, string INPUT_FILENAME, the name of the input file.
//
//    Output, int *M, the number of spatial dimensions.
//
//    Output, int *N, the number of points.
//
{
    *m = file_column_count(input_filename);

    if (*m <= 0) {
        cerr << "\n";
        cerr << "R8MAT_HEADER_READ - Fatal error!\n";
        cerr << "  FILE_COLUMN_COUNT failed.\n";
        *n = -1;
        return;
    }

    *n = file_row_count(input_filename);

    if (*n <= 0) {
        cerr << "\n";
        cerr << "R8MAT_HEADER_READ - Fatal error!\n";
        cerr << "  FILE_ROW_COUNT failed.\n";
        return;
    }

    return;
}
//****************************************************************************80

double *r8mat_indicator(int m, int n)

//****************************************************************************80
//
//  Purpose:
//
//    R8MAT_INDICATOR sets up an "indicator" R8MAT.
//
//  Discussion:
//
//    An R8MAT is an array of R8's.
//
//    The value of each entry suggests its location, as in:
//
//      11  12  13  14
//      21  22  23  24
//      31  32  33  34
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    25 January 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int M, the number of rows of the matrix.
//    M must be positive.
//
//    Input, int N, the number of columns of the matrix.
//    N must be positive.
//
//    Output, double TABLE[M*N], the indicator matrix.
//
{
    int fac;
    int i;
    int j;
    double *table;

    table = new double[m * n];

    fac = (int) pow(10.0, (i4_log_10(n) + 1));

    for (i = 1; i <= m; i++) {
        for (j = 1; j <= n; j++) {
            table[i - 1 + (j - 1) * m] = (double) (fac * i + j);
        }
    }

    return table;
}
//****************************************************************************80

void r8mat_print(int m, int n, double a[], string title)

//****************************************************************************80
//
//  Purpose:
//
//    R8MAT_PRINT prints an R8MAT, with an optional title.
//
//  Discussion:
//
//    An R8MAT is an array of R8's.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    29 August 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int M, the number of rows in A.
//
//    Input, int N, the number of columns in A.
//
//    Input, double A[M*N], the M by N matrix.
//
//    Input, string TITLE, a title to be printed.
//
{
    r8mat_print_some(m, n, a, 1, 1, m, n, title);

    return;
}
//****************************************************************************80

void r8mat_print_some(int m, int n, double a[], int ilo, int jlo, int ihi,
        int jhi, string title)

//****************************************************************************80
//
//  Purpose:
//
//    R8MAT_PRINT_SOME prints some of an R8MAT.
//
//  Discussion:
//
//    An R8MAT is an array of R8's.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    09 April 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int M, the number of rows of the matrix.
//    M must be positive.
//
//    Input, int N, the number of columns of the matrix.
//    N must be positive.
//
//    Input, double A[M*N], the matrix.
//
//    Input, int ILO, JLO, IHI, JHI, designate the first row and
//    column, and the last row and column to be printed.
//
//    Input, string TITLE, a title for the matrix.
{
#define INCX 5

    int i;
    int i2hi;
    int i2lo;
    int j;
    int j2hi;
    int j2lo;

    if (0 < s_len_trim(title)) {
        cout << "\n";
        cout << title << "\n";
    }
    //
    //  Print the columns of the matrix, in strips of 5.
    //
    for (j2lo = jlo; j2lo <= jhi; j2lo = j2lo + INCX) {
        j2hi = j2lo + INCX - 1;
        j2hi = i4_min(j2hi, n);
        j2hi = i4_min(j2hi, jhi);

        cout << "\n";
        //
        //  For each column J in the current range...
        //
        //  Write the header.
        //
        cout << "  Col:    ";
        for (j = j2lo; j <= j2hi; j++) {
            cout << setw(7) << j << "       ";
        }
        cout << "\n";
        cout << "  Row\n";
        cout << "  ---\n";
        //
        //  Determine the range of the rows in this strip.
        //
        i2lo = i4_max(ilo, 1);
        i2hi = i4_min(ihi, m);

        for (i = i2lo; i <= i2hi; i++) {
            //
            //  Print out (up to) 5 entries in row I, that lie in the current strip.
            //
            cout << setw(5) << i << "  ";
            for (j = j2lo; j <= j2hi; j++) {
                cout << setw(12) << a[i - 1 + (j - 1) * m] << "  ";
            }
            cout << "\n";
        }

    }

    return;
#undef INCX
}
//****************************************************************************80

double *r8mat_read(string input_filename, int *m, int *n)

//****************************************************************************80
//
//  Purpose:
//
//    R8MAT_READ reads information from an R8MAT file.
//
//  Discussion:
//
//    An R8MAT is an array of R8's.
//
//    The file is assumed to contain one record per line.
//
//    Records beginning with '#' are comments, and are ignored.
//    Blank lines are also ignored.
//
//    Each line that is not ignored is assumed to contain exactly (or at least)
//    M real numbers, representing the coordinates of a point.
//
//    There are assumed to be exactly (or at least) N such records.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 February 2009
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, string INPUT_FILENAME, the name of the input file.
//
//    Output, int *M, the number of spatial dimensions.
//
//    Output, int *N, the number of points.  The program
//    will stop reading data once N values have been read.
//
//    Output, double R8MAT_READ[M*N], the table data.
//
{
    double *table;

    r8mat_header_read(input_filename, m, n);

    table = r8mat_data_read(input_filename, *m, *n);

    return table;
}
//****************************************************************************80

void r8mat_transpose_print(int m, int n, double a[], string title)

//****************************************************************************80
//
//  Purpose:
//
//    R8MAT_TRANSPOSE_PRINT prints an R8MAT, transposed.
//
//  Discussion:
//
//    An R8MAT is an array of R8's.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    11 August 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int M, N, the number of rows and columns.
//
//    Input, double A[M*N], an M by N matrix to be printed.
//
//    Input, string TITLE, an optional title.
//
{
    r8mat_transpose_print_some(m, n, a, 1, 1, m, n, title);

    return;
}
//****************************************************************************80

void r8mat_transpose_print_some(int m, int n, double a[], int ilo, int jlo,
        int ihi, int jhi, string title)

//****************************************************************************80
//
//  Purpose:
//
//    R8MAT_TRANSPOSE_PRINT_SOME prints some of an R8MAT, transposed.
//
//  Discussion:
//
//    An R8MAT is an array of R8's.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    11 August 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int M, N, the number of rows and columns.
//
//    Input, double A[M*N], an M by N matrix to be printed.
//
//    Input, int ILO, JLO, the first row and column to print.
//
//    Input, int IHI, JHI, the last row and column to print.
//
//    Input, string TITLE, an optional title.
//
{
#define INCX 5

    int i;
    int i2;
    int i2hi;
    int i2lo;
    int inc;
    int j;
    int j2hi;
    int j2lo;

    if (0 < s_len_trim(title)) {
        cout << "\n";
        cout << title << "\n";
    }

    for (i2lo = i4_max(ilo, 1); i2lo <= i4_min(ihi, m); i2lo = i2lo + INCX) {
        i2hi = i2lo + INCX - 1;
        i2hi = i4_min(i2hi, m);
        i2hi = i4_min(i2hi, ihi);

        inc = i2hi + 1 - i2lo;

        cout << "\n";
        cout << "  Row: ";
        for (i = i2lo; i <= i2hi; i++) {
            cout << setw(7) << i << "       ";
        }
        cout << "\n";
        cout << "  Col\n";

        j2lo = i4_max(jlo, 1);
        j2hi = i4_min(jhi, n);

        for (j = j2lo; j <= j2hi; j++) {
            cout << setw(5) << j << " ";
            for (i2 = 1; i2 <= inc; i2++) {
                i = i2lo - 1 + i2;
                cout << setw(14) << a[(i - 1)+(j - 1) * m];
            }
            cout << "\n";
        }
    }

    return;
#undef INCX
}
//****************************************************************************80

double *r8mat_uniform_01(int m, int n, int *seed)

//****************************************************************************80
//
//  Purpose:
//
//    R8MAT_UNIFORM_01 returns a unit pseudorandom R8MAT.
//
//  Discussion:
//
//    An R8MAT is an array of R8's.
//
//    This routine implements the recursion
//
//      seed = 16807 * seed mod ( 2**31 - 1 )
//      unif = seed / ( 2**31 - 1 )
//
//    The integer arithmetic never requires more than 32 bits,
//    including a sign bit.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    03 October 2005
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Paul Bratley, Bennett Fox, Linus Schrage,
//    A Guide to Simulation,
//    Springer Verlag, pages 201-202, 1983.
//
//    Bennett Fox,
//    Algorithm 647:
//    Implementation and Relative Efficiency of Quasirandom
//    Sequence Generators,
//    ACM Transactions on Mathematical Software,
//    Volume 12, Number 4, pages 362-376, 1986.
//
//    Peter Lewis, Allen Goodman, James Miller,
//    A Pseudo-Random Number Generator for the System/360,
//    IBM Systems Journal,
//    Volume 8, pages 136-143, 1969.
//
//  Parameters:
//
//    Input, int M, N, the number of rows and columns.
//
//    Input/output, int *SEED, the "seed" value.  Normally, this
//    value should not be 0.  On output, SEED has 
//    been updated.
//
//    Output, double R8MAT_UNIFORM_01[M*N], a matrix of pseudorandom values.
//
{
    int i;
    int j;
    int k;
    double *r;

    if (*seed == 0) {
        cerr << "\n";
        cerr << "R8MAT_UNIFORM_01 - Fatal error!\n";
        cerr << "  Input value of SEED = 0.\n";
        exit(1);
    }

    r = new double[m * n];

    for (j = 0; j < n; j++) {
        for (i = 0; i < m; i++) {
            k = *seed / 127773;

            *seed = 16807 * (*seed - k * 127773) - k * 2836;

            if (*seed < 0) {
                *seed = *seed + 2147483647;
            }

            r[i + j * m] = (double) (*seed) * 4.656612875E-10;
        }
    }

    return r;
}
//****************************************************************************80

void r8mat_write(string output_filename, int m, int n, double table[])

//****************************************************************************80
//
//  Purpose:
//
//    R8MAT_WRITE writes an R8MAT file.
//
//  Discussion:
//
//    An R8MAT is an array of R8's.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    29 June 2009
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, string OUTPUT_FILENAME, the output filename.
//
//    Input, int M, the spatial dimension.
//
//    Input, int N, the number of points.
//
//    Input, double TABLE[M*N], the table data.
//
{
    int i;
    int j;
    ofstream output;
    //
    //  Open the file.
    //
    output.open(output_filename.c_str());

    if (!output) {
        cerr << "\n";
        cerr << "R8MAT_WRITE - Fatal error!\n";
        cerr << "  Could not open the output file.\n";
        return;
    }
    //
    //  Write the data.
    //
    for (j = 0; j < n; j++) {
        for (i = 0; i < m; i++) {
            output << "  " << setw(24) << setprecision(16) << table[i + j * m];
        }
        output << "\n";
    }
    //
    //  Close the file.
    //
    output.close();

    return;
}
//****************************************************************************80

int s_len_trim(string s)

//****************************************************************************80
//
//  Purpose:
//
//    S_LEN_TRIM returns the length of a string to the last nonblank.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    05 July 2009
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, string S, a string.
//
//    Output, int S_LEN_TRIM, the length of the string to the last nonblank.
//    If S_LEN_TRIM is 0, then the string is entirely blank.
//
{
    int n;

    n = s.length();

    while (0 < n) {
        if (s[n - 1] != ' ') {
            return n;
        }
        n = n - 1;
    }

    return n;
}
//****************************************************************************80

int s_to_i4(string s, int *last, bool *error)

//****************************************************************************80
//
//  Purpose:
//
//    S_TO_I4 reads an I4 from a string.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    05 July 2009
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, string S, a string to be examined.
//
//    Output, int *LAST, the last character of S used to make IVAL.
//
//    Output, bool *ERROR is TRUE if an error occurred.
//
//    Output, int *S_TO_I4, the integer value read from the string.
//    If the string is blank, then IVAL will be returned 0.
//
{
    char c;
    int i;
    int isgn;
    int istate;
    int ival;

    *error = false;
    istate = 0;
    isgn = 1;
    i = 0;
    ival = 0;

    for (;;) {
        c = s[i];
        i = i + 1;
        //
        //  Haven't read anything.
        //
        if (istate == 0) {
            if (c == ' ') {
            } else if (c == '-') {
                istate = 1;
                isgn = -1;
            } else if (c == '+') {
                istate = 1;
                isgn = +1;
            } else if ('0' <= c && c <= '9') {
                istate = 2;
                ival = c - '0';
            } else {
                *error = true;
                return ival;
            }
        }            //
            //  Have read the sign, expecting digits.
            //
        else if (istate == 1) {
            if (c == ' ') {
            } else if ('0' <= c && c <= '9') {
                istate = 2;
                ival = c - '0';
            } else {
                *error = true;
                return ival;
            }
        }            //
            //  Have read at least one digit, expecting more.
            //
        else if (istate == 2) {
            if ('0' <= c && c <= '9') {
                ival = 10 * (ival) + c - '0';
            } else {
                ival = isgn * ival;
                *last = i - 1;
                return ival;
            }

        }
    }
    //
    //  If we read all the characters in the string, see if we're OK.
    //
    if (istate == 2) {
        ival = isgn * ival;
        *last = s_len_trim(s);
    } else {
        *error = true;
        *last = 0;
    }

    return ival;
}
//****************************************************************************80

bool s_to_i4vec(string s, int n, int ivec[])

//****************************************************************************80
//
//  Purpose:
//
//    S_TO_I4VEC reads an I4VEC from a string.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    05 July 2009
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, string S, the string to be read.
//
//    Input, int N, the number of values expected.
//
//    Output, int IVEC[N], the values read from the string.
//
//    Output, bool S_TO_I4VEC, is TRUE if an error occurred.
//
{
    int begin;
    bool error;
    int i;
    int lchar;
    int length;

    begin = 0;
    length = s.length();
    error = 0;

    for (i = 0; i < n; i++) {
        ivec[i] = s_to_i4(s.substr(begin, length), &lchar, &error);

        if (error) {
            return error;
        }
        begin = begin + lchar;
        length = length - lchar;
    }

    return error;
}
//****************************************************************************80

double s_to_r8(string s, int *lchar, bool *error)

//****************************************************************************80
//
//  Purpose:
//
//    S_TO_R8 reads an R8 from a string.
//
//  Discussion:
//
//    This routine will read as many characters as possible until it reaches
//    the end of the string, or encounters a character which cannot be
//    part of the real number.
//
//    Legal input is:
//
//       1 blanks,
//       2 '+' or '-' sign,
//       2.5 spaces
//       3 integer part,
//       4 decimal point,
//       5 fraction part,
//       6 'E' or 'e' or 'D' or 'd', exponent marker,
//       7 exponent sign,
//       8 exponent integer part,
//       9 exponent decimal point,
//      10 exponent fraction part,
//      11 blanks,
//      12 final comma or semicolon.
//
//    with most quantities optional.
//
//  Example:
//
//    S                 R
//
//    '1'               1.0
//    '     1   '       1.0
//    '1A'              1.0
//    '12,34,56'        12.0
//    '  34 7'          34.0
//    '-1E2ABCD'        -100.0
//    '-1X2ABCD'        -1.0
//    ' 2E-1'           0.2
//    '23.45'           23.45
//    '-4.2E+2'         -420.0
//    '17d2'            1700.0
//    '-14e-2'         -0.14
//    'e2'              100.0
//    '-12.73e-9.23'   -12.73 * 10.0**(-9.23)
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    05 July 2009
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, string S, the string containing the
//    data to be read.  Reading will begin at position 1 and
//    terminate at the end of the string, or when no more
//    characters can be read to form a legal real.  Blanks,
//    commas, or other nonnumeric data will, in particular,
//    cause the conversion to halt.
//
//    Output, int *LCHAR, the number of characters read from
//    the string to form the number, including any terminating
//    characters such as a trailing comma or blanks.
//
//    Output, bool *ERROR, is true if an error occurred.
//
//    Output, double S_TO_R8, the real value that was read from the string.
//
{
    char c;
    int ihave;
    int isgn;
    int iterm;
    int jbot;
    int jsgn;
    int jtop;
    int nchar;
    int ndig;
    double r;
    double rbot;
    double rexp;
    double rtop;
    char TAB = 9;

    nchar = s_len_trim(s);
    *error = false;
    r = 0.0;
    *lchar = -1;
    isgn = 1;
    rtop = 0.0;
    rbot = 1.0;
    jsgn = 1;
    jtop = 0;
    jbot = 1;
    ihave = 1;
    iterm = 0;

    for (;;) {
        c = s[*lchar + 1];
        *lchar = *lchar + 1;
        //
        //  Blank or TAB character.
        //
        if (c == ' ' || c == TAB) {
            if (ihave == 2) {
            } else if (ihave == 6 || ihave == 7) {
                iterm = 1;
            } else if (1 < ihave) {
                ihave = 11;
            }
        }            //
            //  Comma.
            //
        else if (c == ',' || c == ';') {
            if (ihave != 1) {
                iterm = 1;
                ihave = 12;
                *lchar = *lchar + 1;
            }
        }            //
            //  Minus sign.
            //
        else if (c == '-') {
            if (ihave == 1) {
                ihave = 2;
                isgn = -1;
            } else if (ihave == 6) {
                ihave = 7;
                jsgn = -1;
            } else {
                iterm = 1;
            }
        }            //
            //  Plus sign.
            //
        else if (c == '+') {
            if (ihave == 1) {
                ihave = 2;
            } else if (ihave == 6) {
                ihave = 7;
            } else {
                iterm = 1;
            }
        }            //
            //  Decimal point.
            //
        else if (c == '.') {
            if (ihave < 4) {
                ihave = 4;
            } else if (6 <= ihave && ihave <= 8) {
                ihave = 9;
            } else {
                iterm = 1;
            }
        }            //
            //  Exponent marker.
            //
        else if (ch_eqi(c, 'E') || ch_eqi(c, 'D')) {
            if (ihave < 6) {
                ihave = 6;
            } else {
                iterm = 1;
            }
        }            //
            //  Digit.
            //
        else if (ihave < 11 && '0' <= c && c <= '9') {
            if (ihave <= 2) {
                ihave = 3;
            } else if (ihave == 4) {
                ihave = 5;
            } else if (ihave == 6 || ihave == 7) {
                ihave = 8;
            } else if (ihave == 9) {
                ihave = 10;
            }

            ndig = ch_to_digit(c);

            if (ihave == 3) {
                rtop = 10.0 * rtop + (double) ndig;
            } else if (ihave == 5) {
                rtop = 10.0 * rtop + (double) ndig;
                rbot = 10.0 * rbot;
            } else if (ihave == 8) {
                jtop = 10 * jtop + ndig;
            } else if (ihave == 10) {
                jtop = 10 * jtop + ndig;
                jbot = 10 * jbot;
            }

        }            //
            //  Anything else is regarded as a terminator.
            //
        else {
            iterm = 1;
        }
        //
        //  If we haven't seen a terminator, and we haven't examined the
        //  entire string, go get the next character.
        //
        if (iterm == 1 || nchar <= *lchar + 1) {
            break;
        }

    }
    //
    //  If we haven't seen a terminator, and we have examined the
    //  entire string, then we're done, and LCHAR is equal to NCHAR.
    //
    if (iterm != 1 && (*lchar) + 1 == nchar) {
        *lchar = nchar;
    }
    //
    //  Number seems to have terminated.  Have we got a legal number?
    //  Not if we terminated in states 1, 2, 6 or 7!
    //
    if (ihave == 1 || ihave == 2 || ihave == 6 || ihave == 7) {
        *error = true;
        return r;
    }
    //
    //  Number seems OK.  Form it.
    //
    if (jtop == 0) {
        rexp = 1.0;
    } else {
        if (jbot == 1) {
            rexp = pow(10.0, jsgn * jtop);
        } else {
            rexp = jsgn * jtop;
            rexp = rexp / jbot;
            rexp = pow(10.0, rexp);
        }

    }

    r = isgn * rexp * rtop / rbot;

    return r;
}
//****************************************************************************80

bool s_to_r8vec(string s, int n, double rvec[])

//****************************************************************************80
//
//  Purpose:
//
//    S_TO_R8VEC reads an R8VEC from a string.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    05 July 2009
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, string S, the string to be read.
//
//    Input, int N, the number of values expected.
//
//    Output, double RVEC[N], the values read from the string.
//
//    Output, bool S_TO_R8VEC, is true if an error occurred.
//
{
    int begin;
    bool error;
    int i;
    int lchar;
    int length;

    begin = 0;
    length = s.length();
    error = 0;

    for (i = 0; i < n; i++) {
        rvec[i] = s_to_r8(s.substr(begin, length), &lchar, &error);

        if (error) {
            return error;
        }
        begin = begin + lchar;
        length = length - lchar;
    }

    return error;
}
//****************************************************************************80

int s_word_count(string s)

//****************************************************************************80
//
//  Purpose:
//
//    S_WORD_COUNT counts the number of "words" in a string.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    05 July 2009
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, string S, the string to be examined.
//
//    Output, int S_WORD_COUNT, the number of "words" in the string.
//    Words are presumed to be separated by one or more blanks.
//
{
    bool blank;
    int char_count;
    int i;
    int word_count;

    word_count = 0;
    blank = true;

    char_count = s.length();

    for (i = 0; i < char_count; i++) {
        if (isspace(s[i]) || s[i] == ',') {
            blank = true;
        } else if (blank) {
            word_count = word_count + 1;
            blank = false;
        }
    }

    return word_count;
}
//****************************************************************************80

void timestamp()

//****************************************************************************80
//
//  Purpose:
//
//    TIMESTAMP prints the current YMDHMS date as a time stamp.
//
//  Example:
//
//    31 May 2001 09:45:54 AM
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    08 July 2009
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    None
//
{
#define TIME_SIZE 40

    static char time_buffer[TIME_SIZE];
    const struct std::tm *tm_ptr;
    std::time_t now;

    now = std::time(NULL);
    tm_ptr = std::localtime(&now);

    std::strftime(time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm_ptr);

    std::cout << time_buffer << "\n";

    return;
#undef TIME_SIZE
}

extern "C" int table_column_count(char *input_filename) {
    return (file_column_count(input_filename));
}

extern "C" int table_row_count(char *input_filename) {
    return (file_row_count(input_filename));
}

extern "C" int *table_read_i4(char *input_filename, int m, int n) {
    return (i4mat_data_read(input_filename, m, n));
}

extern "C" double *table_read_r8(char *input_filename, int m, int n) {
    return (r8mat_data_read(input_filename, m, n));
}

extern "C" float *table_read_r4(char *input_filename, int m, int n) {
    return (r4mat_data_read(input_filename, m, n));
}

extern "C" void table_free_i4(int *table) {
    delete [] table;
}

extern "C" void table_free_r8(double *table) {
    delete [] table;
}

extern "C" void table_free_r4(float *table) {
    delete [] table;
}

